function [Output_struc,inputs_struc] = solve_counterfactual(inputs_struc,fns_struc)

% Takes the supply of skill and target pcGDP and solves for the equilibrium
% in the counterfactual adjusting theta, relative productivies and average
% productivity to match these targets while maintaining the w_S and
% relative prices.

w_U = inputs_struc.w_U; w_S_target = inputs_struc.w_S_target;sigma_us = inputs_struc.sigma_us;
prod_grid_size=inputs_struc.prod_grid_size; no_q = inputs_struc.no_q; A_S = inputs_struc.A_S; 
L_U_by_L_S_q_target = inputs_struc.L_U_by_L_S_q_target; tol = inputs_struc.tol;
A_q_i = inputs_struc.A_q_i; price_size_slope_target = inputs_struc.price_size_slope_target;
size_q_target = inputs_struc.size_q_target; size_hist_cutoffs_graph = inputs_struc.size_hist_cutoffs_graph;
p_q_i_calib = inputs_struc.p_q_i_calib; pdf_A_q_i = inputs_struc.pdf_A_q_i;
N = inputs_struc.N; pcgdp_calib_own_prices = inputs_struc.pcgdp_calib_own_prices;

MC_L_part_fn = fns_struc.MC_L_part_fn;

prod_mult_guess = inputs_struc.pcgdp_target;

% Bisection over thetas to meet wage premia
error_theta = 1;
iter_theta = 0;
min_theta = 0;
max_theta = 1;

tic
while error_theta>tol*10 && iter_theta<30 

    % Bisection over theta of lowest quality with all higher qualities
    % chosen to to match skilled to unskilled ratio across qualities
    theta_bisec = (min_theta + max_theta)/2;
    theta_bisec_vec = 1./(1 + ((L_U_by_L_S_q_target(1,1)./L_U_by_L_S_q_target(2:end,1)).^(1/sigma_us)).*((1-theta_bisec)/theta_bisec));                  % Share parameter of unskilled in high quality
    theta_bisec_vec = [theta_bisec;theta_bisec_vec];
    iter_theta = iter_theta + 1;
    theta_q_i_bisec = repmat(theta_bisec_vec,[1,prod_grid_size]);
    inputs_struc.theta_q_i = theta_q_i_bisec;

    % The entry denomination given thetas and target w_S
    S_to_U_ratio_q_bisec = ((w_U/w_S_target*(1-theta_bisec_vec)./theta_bisec_vec).^(sigma_us)).*(A_S.^(sigma_us-1));
    entry_sh_S_q_bisec = 1 - 1./(S_to_U_ratio_q_bisec + 1);
    inputs_struc.entry_sh_S_q = entry_sh_S_q_bisec;

    % Given thetas, figure out the prod_mult to match pcGDP. Inside the
    % function, i make sure that relative priced dont change
    prod_mult_guess = C_theta_prodmult_fn(prod_mult_guess,inputs_struc,fns_struc); 

    % Just solve for equi with the calibrated prod_mult
    price_ratio_target_q  = exp(price_size_slope_target*log(size_q_target));  
    price_ratio_target_q = price_ratio_target_q/price_ratio_target_q(1,1);      % Target for aggregate prices across qualities.
    MC_L_part = MC_L_part_fn(w_U,w_S_target,theta_bisec_vec,sigma_us,A_S);
    A_scale_q_bisec = (MC_L_part(1,1)./MC_L_part)./price_ratio_target_q; %ones(no_q,1); %
    A_q_i_bisec = repmat(A_q_i(1,:),[no_q,1]).*repmat(A_scale_q_bisec,1,prod_grid_size);
    inputs_struc.A_q_i = A_q_i_bisec*prod_mult_guess ;
    Output_struc = solve_equi(inputs_struc,fns_struc);
    inputs_struc.A_q_i = A_q_i;


    w_S_ch_theta_c = Output_struc.w_S ;
    error_theta = abs(w_S_ch_theta_c/w_S_target - 1);

    if w_S_ch_theta_c>w_S_target
        min_theta = theta_bisec;
    else
        max_theta = theta_bisec;
    end

end

toc

% Computing Size Distribution
l_U_q_i = Output_struc.l_U_q_i;
l_S_q_i = Output_struc.l_S_q_i;
size_q_i_C = l_U_q_i + l_S_q_i;
mass_firms_q_i_C = inputs_struc.pdf_A_q_i.*repmat(Output_struc.M_q,[1,prod_grid_size]);
for j = 1:size(size_hist_cutoffs_graph,1)-1

    indicator_model = (size_q_i_C<=size_hist_cutoffs_graph(j));
    size_hist_cdf_model_C_across_St(1,j)            = sum(sum(indicator_model.*size_q_i_C.*mass_firms_q_i_C))/ sum(sum(size_q_i_C.*mass_firms_q_i_C));

end

size_hist_cdf_model_C_across_St(1,j+1) = 1;

inputs_struc.A_scale_q_bisec = A_scale_q_bisec;
inputs_struc.prod_mult_guess = prod_mult_guess;

d_q_i_C=Output_struc.d_q_i;
M_q_C=Output_struc.M_q;
pcgdp_prodvec_C_calib_prices = sum(sum(p_q_i_calib.*d_q_i_C.*pdf_A_q_i.*repmat(M_q_C,[1,size(p_q_i_calib,2)])))/N;
pcgdp_calib_rel = pcgdp_prodvec_C_calib_prices/pcgdp_calib_own_prices;

% Stuff which will be written into excel
write_excel = [inputs_struc.min_q;inputs_struc.h;inputs_struc.prod_mult_guess;pcgdp_calib_rel;Output_struc.w_S;Output_struc.L_U_clearing;Output_struc.L_S_clearing];
rel_price_write = Output_struc.P_q_norm;
rel_price_write = rel_price_write./rel_price_write(1);
avg_size_write = (Output_struc.L_U_q + Output_struc.L_S_q)./Output_struc.M_q;
theta_write = inputs_struc.theta_q_i;
theta_write = theta_write(:,1);
rel_prod_write = inputs_struc.A_q_i;
rel_prod_write = rel_prod_write(:,1)./rel_prod_write(1,1);
sh_S_entry_write = inputs_struc.entry_sh_S_q;
write_excel = [write_excel;0;0;rel_price_write; avg_size_write; theta_write; rel_prod_write ;sh_S_entry_write];

Output_struc.pcgdp_calib_rel = pcgdp_calib_rel;
Output_struc.write_excel = write_excel;
Output_struc.size_hist_cdf_model_C = size_hist_cdf_model_C_across_St;



end
